Advanced stability analysis of a fractional delay differential system with stochastic phenomena using spectral collocation method

In recent years, there has been a growing interest in incorporating fractional calculus into stochastic delay systems due to its ability to model complex phenomena with uncertainties and memory effects. The fractional stochastic delay differential equations are conventional in modeling such complex dynamical systems around various applied fields. The present study addresses a novel spectral approach to demonstrate the stability behavior and numerical solution of the systems characterized by stochasticity along with fractional derivatives and time delay. By bridging the gap between fractional calculus, stochastic processes, and spectral analysis, this work contributes to the field of fractional dynamics and enriches the toolbox of analytical tools available for investigating the stability of systems with delays and uncertainties. To illustrate the practical implications and validate the theoretical findings of our approach, some numerical simulations are presented.

Fractional derivatives ( D α and D β ) D α U(t) denotes the fractional derivative of U(t) of a fractional order α , that capturing the long-range dependen- cies and memory effects.Where D β V (t) denotes the fractional derivative of V(t) of a fractional order β.

Deterministic components
To represent the deterministic drift terms along with time delay for U(t) and V(t) are a(t)U(t) + b(t)V (t − τ ) and c(t)U(t) + d(t)V (t − τ ) respectively.These terms require the deterministic expansion of the proposed variables.

Stochastic components
The stochastic components of the system of equations are represented by σ 1 (t)dW(t) and σ 2 (t)dW(t) .These terms present the random fluctuations into the dynamics of U(t) and V(t) respectively.The proposed FSDDEs system have been models the dynamic interaction between both variables U(t) and V(t).However, their evolution is convinced by both stochastic noise represented by the Wiener process dW(t) and fractional derivatives capturing memory effects.All the coefficients a(t), b(t), c(t) and d(t) represent the nature and strength of the interaction between the deterministic evolution and their variables.
The study of the proposed such systems is very significant in various scientific research fields, such as physics, finance, and biology, where random fluctuations and memory effects play a significant role in conforming the observed dynamics.

Formulation with the Caputo fractional derivative
The Caputo derivative applied to the proposed FSDDEs system introduced earlier section.The Caputo fractional derivative (CFD) is a one of the common choice when modeling physical systems with fractional dynamics, as it better to apply for the well-posed initial value problems.
In general, the Caputo fractional derivative of a function f(t) with order α is defined by: where D C α denotes the Caputo type derivative, α is the fractional order of the given derivative, where 0 < α ≤ 1 , Ŵ denotes the gamma function , f ′ (s) denotes the first derivative of f .Now, rewrite the system of two Caputo FSDDEs as follows: (1) • CFD for U(t): D C α U(t) denotes the CFD of U(t) with fractional order α .The proposed fractional derivative captures the memory effect in the evolution of the variable U(t), allowing it to the rate of change of U(t) and depend on past values.
• CFD for V(t): D C β V (t) denotes the CFD of V(t) of order β .Similarly, the above fractional derivative also captures the memory effect in the evolution of the variable V(t).
The presence of CFD in the system implies that the variables U(t) and V(t) dynamics exhibition the fractional order memory effects, which make the system capable for the systems modeling with anomalous diffusion and long-range dependencies.
The proposed system also includes stochastic components, interaction coefficients and deterministic drift terms, as explained in the previous section.The choice of CFD guarantees that system is suitable and well-posed for solving the initial value problems, for modeling the complex systems with stochastic influences and fractional dynamics it making a powerful tool.

Stability analysis of FSDDEs
The stability analysis of FSDDEs using the spectral collocation technique concerns interrogating the behavior of the proposed system eigenvalues to determine whether it offerings stable or an unstable dynamics.This approach along with principles of fractional calculus with spectral theory to assess the system stability.Here, in present research work we provide a step-by-step mathematical approach to managing stability analysis by using the spectral collocation method.

Linearization of FSDDEs
Linearization of FSDDEs is a difficult step in stability analysis.Linearization allows us to simplify complex nonlinear system around to analyze their behavior and an equilibrium point.In this detailed discussion, we'll explore the process of linearization for given FSDDEs and understand the mathematical philosophy behind it.

Starting point: nonlinear FSDDEs
We begin with a system of nonlinear FSDDEs in the general form: System Eq. ( 4) involve the fractional derivatives D α and D β , where f and g are nonlinear functions and sto- chastic terms along with Wiener processes σ 1 (t)dW(t) and σ 2 (t)dW(t).

Equilibrium point
To analyze stability analysis, we first to identify an equilibrium point (U * , V * ) , for this take both the time deriva- tives of D α U(t) and D β V (t) are zero: Note that, these equilibrium values represent a stable state for the proposed system.

Linearization process
Around the equilibrium point, we first introduce small perturbations: Substituting the perturbations Eq. ( 5) into the FSDDEs given in Eq. ( 4) and using Taylor series expansions, we keep only the linear terms in the perturbations: When the perturbations are small, the higher-order terms like D α U(t) 2 and D β V (t) 2 can be ignored, because they are very small as compared to linear terms.

Linear system
After successfully linearized the proposed system FSDDEs given in Eq. ( 4), the resulting equations are a system of linear stochastic delay differential equations (LSDDEs) that characterize the behavior of small perturbations U(t) and V(t) around the present equilibrium point.

The Fourier transform in the analysis of FSDDEs
Fourier transform is one of the powerful mathematical tool used in the analysis of system of FSDDEs.Fourier transformation empowers us to work in frequency domain, where the analysis of FSDDEs and linearization become more controllable.We'll explore how to apply the Fourier transform to a system of FSDDEs and its inferences for stability analysis.

Fourier transform basics
For any function f(t), the continuous Fourier transform f (ω) is given as: here ω is the angular frequency and f (ω) is the complex-valued frequency domain that represent f(t).

Applying the Fourier transform to FSDDEs
We'll continue with the linearized equilibrium points given in Eq. ( 5) and apply the Fourier Transform.We'll use F to denote the Fourier Transform: where U * (ω) and V * (ω) denote the respective Fourier transforms of U(t) and V(t).Similarly, the variables u(ω) and v(ω) are Fourier transforms of the perturbations U(t) and V(t).
Next, we apply the Fourier Transform to the linearized FSDDEs Eq. (1).Given that a, b, c, and d are constants, we can represent linearized system of FSDDEs in the frequency domain:

Fourier transform of derivatives
To manipulate both the fractional derivatives D α and D β in frequency domain, we must use the Fourier transfor- mation properties.Specifically, if f (ω) is the Fourier transform of a function f(t), then the fractional derivative D α f (t) is given by: the operator F denotes the Fourier transformation.Now using the above property to the system of FSDDEs given in Eq. ( 8) yields: where Û(ω) , V (ω) , and Ŵ(ω) are the Fourier transforms of the functions U(t), V(t), and W(t), respectively.

Eigenvalue analysis in the frequency domain
In frequency domain of the Fourier-transformed FSDDEs, the eigenvalue analysis involves the finding the eigenvalues of the proposed system matrix to impose stability of the linearized FSDDEs Eq. ( 9) in matrix form: The above system is a linear system in a frequency domain, where the functions Û(ω) and V (ω) are both unknowns, where the matrix on a left hand side denotes the system matrix A(ω).Characteristic equation Now we are find eigenvalues of the proposed system matrix A(ω) , for this to calculate the given characteristic equation: here I is identity matrix.In this case, the given characteristic equation having order 2x2 and define as:

Eigenvalues and stability
Here, to decide the stability of the proposed system, we must need to determine the eigenvalues of system matrix in a frequency domain.In this case, the eigenvalues corresponding to the roots of the characteristic equation which derived by from the above linearized FSDDEs in the frequency domain: To assessing the stability analysis, we interrogate the given real parts of the eigenvalues, that is Re( ) .Then the stability benchmarks are as follows: If all the eigenvalues are negative therefore, Re( ) < 0 , then the proposed system is stable asymptotically.Basically, it shows that all perturbations decay due to time, and the system goes to its equilibrium point.
Conversely, for any eigenvalue, if Re( ) > 0 , then the system is unstable.In the above case, the system has at least one of the perturbation grows over time, indicating divergence from the equilibrium point and instability.
At the same time, if there are both negative and positive real parts through the eigenvalues, in this situation the system shows the complex dynamics.Some perturbations decay, while others grow, that's leading to complex behavior.
To find the eigenvalues , we need to solve the proposed characteristic equation symbolically or numerically for frequency ranges and specific parameter values.The above eigenvalues depend on the values of a, b, c, d, α , β , and ω.
However, keep in mind that the stability of proposed system can change that depending on frequency of perturbations and the parameter values.Therefore, operating stability analysis over frequencies and using different parameter ranges which can provide a comprehensive understanding of a system's behavior.

Inverse Fourier transform
Now to find the inverse Fourier transform, in order to return of the linearized FSDDEs to time domain, for this we need to apply the inverse Fourier transform operator F −1 to the Eq. ( 9) in the frequency domain given in form: Expressions for U(t) and V(t) in the time domain can be obtained using these equations, which can be used for simulations and further analysis.There are various ways in which these expressions can actually look depending on the values of α and β , along with the functional forms of the coefficients a, b, c, and d, and the stochastic process properties Ŵ(ω).

Spectral method
Using the Legendre spectral collocation method to solve fractional stochastic delay differential equations (FSD-DEs) involves several steps.Our method will be to present key equations along with a step-by-step explanation of the process as opposed to presenting the full calculation as a whole.Solving such complex equations would typically require dedicated numerical software or libraries, and the detailed implementation may vary depending on specific parameters and functions involved.Nevertheless, this outline should give you a clear idea of how to approach the problem.
Let's consider a general system of FSDDEs with fractional derivatives of order α and β: State variables U(t) and V(t) are represented by variable f and variable g, standard deviations σ 1 (t) and σ 2 (t) , and variables ξ 1 (t) and ξ 2 (t) are represented by variable f.
Step 1: discretization of time Create discrete points or nodes in the time domain t n , where n = 0, 1, 2, . . ., N .As a result of these nodes, we will be able to construct the collocation method.The Legendre polynomials, denoted as P n (x) , are defined by the following recursive formula: where n is a non-negative integer, and P 0 (x) and P 1 (x) are the initial Legendre polynomials.Express the state variables U(t) and V(t) as expansions in terms of Legendre polynomials: The expansion coefficients a k and b k are determined by the Legendre polynomials used in the expansion, and the scaled time variable τ is determined by the Legendre polynomials used in the expansion [−1, 1].

Step 3: approximation of derivatives
To approximate both the fractional derivatives D α U(t) and D β V (t) using the proposed method.For example, for D α U(t): here the fractional derivative D α P k (τ ) of the Legendre polynomial P k (τ ).
Step 4: discretization of delay terms Also the approximation the delayed term U(t − τ ) at the collocation points t n as given by: where U n denotes the value of state variable U(t) at the corresponding time t = t n − τ.
As above similarly, for the variable V (t − τ ) , we can approximate it at t n as: again V n denotes the value of the state variable V(t) at time t = t n − τ .Now replace the original delay terms in your FSDDEs with their discretized counterparts.For example, your FSDDEs might now look like this: Here, U n and V n are the values of U(t) and V(t) at the corresponding discretization nodes t = t n − τ .The equations now incorporate the delayed values as they were approximated at specific time points.
Step 5: collocation equations Using the collocation method, evaluate the FSDDEs at the collocation points t n , yielding a set of algebraic equa- tions involving the expansion coefficients a k and b k .
Step 6: stochastic terms Incorporate the stochastic terms σ 1 (t)ξ 1 (t) and σ 2 (t)ξ 2 (t) in the system equations, usually through Monte Carlo simulations or numerical methods designed for stochastic processes.
Step 7: solve for coefficients To solve for the coefficients a k and b k , we need to plug the approximations of U(t) and V(t) into the FSDDEs and apply the Legendre spectral collocation method.Here's the mathematical Approach.
Substitute the Legendre polynomial expansions given in Eq. ( 11) into the FSDDEs given in Eq. ( 10): Apply the Legendre spectral collocation method: Evaluate the left-hand side (LHS) and right-hand side (RHS) of the FSDDEs at the collocation points t n within the domain [−1, 1] .The collocation points are typically chosen as the roots of the Legendre polynomial of degree K in the interval [−1, 1].
For example, if we choose the K-th degree roots t n of Legendre polynomial, then the collocation system Eq.( 12) for both the variables U(t) and V(t) would be: Apply numerical techniques such as matrix inversion, iterative methods, or specialized solvers to solve the system of algebraic equations Eq. ( 13) for coefficients a k and b k .The complexity of the problem, linear or non- linear are depending on solvers.
The numerical solutions for the state variables U(t) and V(t) due to desired time domain can be constructed by connecting these coefficients back into the Legendre polynomial expansions, then we have determined the values of a k and b k .
It depends on the computational resources available and the complexity of the problem to determine what numerical methods and tools are used to solve the system of equations resulting from the process.The determination of the coefficients using spectral collocation method is often sophisticated using software like Matlab and the numerical libraries.
Step 8: post-processing After obtaining the coefficients, we can use them to construct the approximate solutions U(t) and V(t) over the desired time period.
The method of Legendre spectral collocation is used in this outline to solve a system of FSDDEs.According to the complexity of the problem and the availability of computational resources, specific implementation details, such as collocation points and Legendre polynomials, may vary.

Practical examples
The spectral method can be a challenging and computationally intensive way of solving complex systems of fractional stochastic delay differential equations.Our purpose here is to present three concrete numerical examples, each with varying degrees of complexity, and explain the mathematical calculations involved.In order to conserve space, we will focus on the key steps of each example.For the simulation of the below examples we use Matlab-15.

Example 1 Linear FSDDEs with fractional orders
Consider a system of linear FSDDEs with fractional orders: This system has fractional orders α = 0.5 and β = 0.3 ; where the initial values are X(0) = 1 and Y (0) = 0.5 .Also the time delay τ = 1 and includes stochastic terms.

Solution Discretize the time domain
Choose the number of basis functions for the spectral method and discretize the time domain.For simplicity, let's assume we have N + 1 collocation points with t n = n/N for n = 0, 1, . . ., N.
Express X(t) and Y(t) using Legendre polynomials: Express X(t) and Y(t) as expansions in terms of Legendre polynomials: Here, a k and b k are the expansion coefficients, and τ = 2t − 1 maps the interval [−1, 1] to the original time domain.
Apply the spectral collocation method Apply the spectral collocation method by evaluating the FSDDEs at the collocation points: For n = 0, 1, . . ., N , we have: where τ n = 2t n − 1.
Solve for the coefficients a k and b k For U(t), the FSDDE becomes: where τ n = 2t n − 1. Substitute the Legendre polynomial expansion for U(t) into this equation and apply the spectral collocation method by evaluating it at the collocation points t n .You'll obtain a set of algebraic equations like this for each n: Do the same for D 0.3 Y (t) using the FSDDE for D 0.3 Y (t n ): Now, we have a system of algebraic equations, one for each collocation point t n : These equations are linear in a k and b k , but there will be N + 1 such equations, one for each collocation point t n .You can solve this system of equations numerically using methods like matrix inversion or iterative solvers to obtain the coefficients a k and b k for each k in the expansion.
In the Legendre polynomial expansions, we plug these coefficients back into the Legendre polynomial expansions to obtain the approximate solutions X(t) and Y(t) over the desired time domain.
Performing the inverse Legendre transform will yield the original time domain values X(t) and Y(t).Based on the Legendre polynomial basis functions a k and b k , the reverse Legendre transform is essentially a weighted sum.
Inverse Transform for X(t): Inverse Transform for Y(t): We calculate X(t) and Y(t) at any given time points in the original time domain by substituting the appropriate values of t into the above system.Depending on the specific problem and the accuracy required, you may need to truncate the series to a finite number of terms (K terms) to approximate the solutions.Keep in mind that the accuracy of the solutions will depend on the choice of the number of basis functions (K) and the collocation points, as well as the precision of the numerical methods used to solve for the coefficients.Adjusting these parameters can help you balance computational efficiency with solution accuracy.

Example 2 Nonlinear FSDDEs with time-varying coefficients
Consider a system of nonlinear FSDDEs with time-varying coefficients:

Nonlinear effects
The term sin(t)X(t − τ ) + 0.1X(t − τ ) 2 in the equation for X(t) introduces nonlinearity.It implies that the rate of change of X(t) depends not only on its own past values ( X(t − τ ) ) but also on the current time (t) and its square ( X(t − τ ) 2 ).In a chemical reactions the nonlinearities can arise from complex reaction kinetics, such as inhibitory processes or autocatalytic.

Time-varying coefficients
These equations have time-varying coefficients, such as sin(t) and cos(t) .In reaction the rates or other parameters of a given system, there can be changes due to time.However, the chemical reactions may be vary with time due to external or temperature influences changes.

Fractional orders
In the system Eq.( 15), the fractional orders (0.7 and 0.5) shows that the long-range memory effects and dependencies are considered by fractional derivatives.In addition to making the system more capable and realistic for modeling phenomena along with memory effect, the given fractional orders indicates the influence of past states on the current dynamics.

Stochasticity
The stochastic noise is denoted by σ 1 ξ 1 (t) and σ 2 ξ 2 (t) .In chemical reactions the stochasticity or external influ- ences can cause and also contribute to it.As a result of this example, chemical species concentrations (X(t) and Y(t)) evolve over time in a physical system, like a chemical reaction.There is a nonlinear behavior in the system due to complex reaction kinetics and parameters that vary over time, such as reaction rates.Modeling real-world phenomena becomes more complex by adding fractional derivatives and stochastic terms to account for memory effects and random fluctuations.
In Fig. 4, we show the solution trajectory of FSDDEs system given in Eq. ( 15) for both X(t) and Y(t).Similarly, in Fig. 5, we show the solution of X(t) and Y(t) for different fractional parameter values α and β .Also in Fig. 6, we draw the solution of FSDDEs model given in (Eq.15) for the stochastic parameter value σ = 0.8.

Example 3 Fractional diffusion-reaction system
In this example, we'll consider a complex system of FSDDEs modeling a fractional diffusion-reaction system: Here, u(x, t) and v(x, t) represent concentrations of chemical species, α = 0.8 and β = 0.6 are fractional orders, D = 1 , k = 0.1 , σ = 0.4 , and = 0.2 are constants, and dW(t) represents Wiener increments in time.( 16) Certainly, let's discuss the graphs at different time points obtained from the simulation of the Fractional Diffusion-Reaction System using Caputo derivatives.In the provided code, we captured the concentrations of u and v at various time steps.Below is a discussion of these graphs at different t values: t = 0: In Fig. 7, at the initial time step, t = 0 , both u(x, t) and v(x, t) start with the specified initial conditions.u has a distribution in the middle of the spatial domain, while v starts with a different distribution.t = 2: In Fig. 8, as time progresses, u(x, t) diffuses outward from its initial concentration center due to the diffusion term ( D • u xx ), while v(x, t) evolves as a result of the reaction and diffusion terms.v starts to react   with u and adjusts its distribution accordingly.In the above Fig.8, we observe that only the u 2 (x, t) and v 2 (x, t) graphs are shown and u 1 (x, t) and v 1 (x, t) graphs are disappear, because the range of u 1 (x, t) and v 1 (x, t) is 1 see in Fig. 7, while the range of u 2 (x, t) and v 2 (x, t) is almost 10 10 ; therefore the smaller range go down and becomes disappear.t = 3: In Fig. 9, at this time step, the diffusion process continues, and the profiles of both u and v continue to change.The concentration profiles of u and v evolve, reflecting the diffusion and reaction processes in the system.Once again in Fig. 9, we observe that only the u 3 (x, t) and v 3 (x, t) graphs are shown and the remaining u 1 (x, t) , v 1 (x, t) , u 2 (x, t) and v 2 (x, t) graphs are disappear, because the range of u 3 (x, t) and v 3 (x, t) is almost 10 20 ; therefore the graphs having smaller range go down and becomes disappear.t = 4: In Fig. 10, by t = 3 , we can observe further spreading of u(x, t) and the impact of the reaction term ( −k • u • v ) on both u and v .The distribution of v adapts to the presence of u .Again in Fig. 10, we see the range of u 4 (x, t) and v 4 (x, t) is almost 10 28 ; therefore the remaining graphs u 3 (x, t) and v 3 (x, t) , u 2 (x, t) and v 2 (x, t) and u 1 (x, t) and v 1 (x, t) having smaller range therefore they go down and becomes disappear.t = 100: In Fig. 11, at this point, u(x, t) and v(x, t) continue to evolve.The reaction between u and v becomes more prominent, leading to changes in their concentration profiles.We clearly observed from Fig. 11, that if we increase the time step, both the chemical species u and v increase the concentration with each other.The impact of stochastic noise ( σ • randn ) also plays a role in the variability of the profiles.
These discussions highlight the dynamic behavior of the system over time.As time progresses, the diffusion, reaction, and stochastic terms interact to shape the concentrations of u and v .The plots provide insights into www.nature.com/scientificreports/ the spatiotemporal evolution of the chemical species in the system, showing how they spread, react, and adapt to each other and to external noise.Solving such complex systems typically requires dedicated numerical libraries and computational resources, as well as a deep understanding of both fractional calculus and stochastic processes.

Challenges and future directions
Fractional stochastic systems pose challenges in terms of analytical solutions and numerical simulations due to the combination of non-integer derivatives and stochastic components.However, magnifications in each computational techniques, that is Monte Carlo simulations, spectral methods, and for fractional calculus the numerical algorithms have made it possible to simulate and analyze such systems with high accuracy.In summery, the fractional stochastic systems provide a powerful techniques for such complex dynamics modeling that exhibit both stochastic characteristics and fractional calculus.The proposed systems having applications in a continue to be a subject of active research and wide range of fields, offering perception into the real-world phenomena's behavior of that conventional models may struggle to capture.

Conclusion
In this research work, we have precisely introduced a novel mathematical approach for a systems constrained by FSDDEs and their stability analysis using the spectral collocation technique a strong mathematical approach.By combining the power of spectral techniques, stochastic processes, and fractional calculus, this study offers expensive perception into the stability features of complex dynamical systems along with delays.Numerical simulations confirm the efficiency and accuracy of our approach, for researchers making it a valuable, practitioners working and worthy tool with FSDDEs in various scientific research fields.This knowledge presents substantially for real-world phenomena to our understanding, offering a foundation for predicting and controlling systems exhibiting stochastic dynamics and along with fractional derivative.This research gives the way for further applications and exploration in different fields, from finance to biology.In future this research work is extended to the analysis to more complex systems and investigate the additional applications in engineering and science. https://doi.org/10.1038/s41598-024-62851-0www.nature.com/scientificreports/ 1σ 2 ξ 2 (t).
ξ 2 (t).Step 2: Legendre polynomial expansion Adrien-Marie Legendre was a French mathematician who invented Legendre polynomials.Mathematicians and physicists use these polynomials to solve differential equations and perform numerical analysis.Several properties make Legendre polynomials valuable tools in mathematical modeling and approximation, as they are defined on the interval [−1, 1].